Inverse simulated annealing for the determination of amorphous structures 
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We present a new and efficient optimization method to determine the structure of disordered 
systems in agreement with available experimental data. Our approach permits the application of 
accurate electronic structure calculations within the structure optimization. The new technique is 
demonstrated within density functional theory by the calculation of a model of amorphous carbon. 
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Amorphous solids can be produced from almost any 
chemical system and are of great interest due to their 
large variety of technologically important applications. 
In addition to conventional silicate glasses, they are, for 
example used in optical waveguides (oxides), plastics (or- 
ganic polymers), solar cells (semiconductors), biomate- 
rials (amorphous metals), xerography and non- volatile 
memory devices (chalcogenides), to name but a few 
Nevertheless, finding their atomic scale structure is still 
a major challenge in material science 0-Q due to the 
absence of lattice periodicity and long-range order char- 
acteristics of a crystalline solid. Many sophisticated mod- 
cling techniques from the field of crystal structure pre- 
diction are based on searching the global minimum in 
the energy landscape for periodic structures (5I— f]~3j] . How- 
ever, an amorphous solid does not correspond to a global, 
but to a local energy minimum, which is energetically 
low enough to stabilize the structure against alternative 
packings and exhibits desirable target properties. 

The most commonly applied computational technique 
to obtain the amorphous structure is to slowly quench 
it from the melt by Monte Carlo (MC)- or Molecular 
Dynamics (MD)-based Simulated Annealing (SA) flU ]. 
However, the lack of exploitable symmetry and, there- 
fore, large number of degrees of freedom, require the 
cooling to be conducted as slowly as possible to deter- 
mine an approximation of the amorphous structure, and 
is, therefore, computationally very demanding. This is 
even more pronounced in conjunction with accurate ab- 
initio electronic structure techni ques , in spite of signifi- 
cant progress in recent years lot 1 161 , allowing for satis- 
factory structure determinations [17H20| . 

Instead of performing an elaborate calculation to ob- 
tain an approximate amorphous model and to assess a 
posteriori how well it matches the experiment, McGreevy 
and coworkers demonstrated that it can be beneficial to 
reverse this procedure, hence the name Reverse Monte 
Carlo (RMC) [2l], [22| ■ Contrary to energy-based mini- 
mization techniques this method aims at directly model- 
ing the structure without invoking any computationally 
expensive potential energy calculation, using only avail- 
able experimental data. Specifically, the available ex- 
perimental data are reproduced simply by minimizing a 



function of the form 
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under variation of the atomic positions R = {r^} us- 
ing the Metropolis Monte Carlo method [23| . In Eq. [I] 
X P (R) and Xp Xp are the calculated and experimental val- 
ues, respectively, of a property p, while w p = l/a^ is a 
weight factor and a p is the experimental uncertainty for 
the corresponding property. 

Even though p can, in principle, be any arbitrary prop- 
erty, in practice, only geometric quantities, obtainable 
from Neutron or X-ray scattering data such as the struc- 
ture factor or the pair correlation function for which 
X P (R) can be evaluated easily and fast, are employed. 
In particular, typically no electronic quantities based 
on accurate electronic structure calculations are utilized, 
which would otherwise be computationally unfeasible. 
While on the one hand, RMC allows for an efficient 
and routine modeling of rather complex disordered struc- 
tures, on the other hand, the resulting models are not 
necessarily physically sensible. It is, therefore, good prac- 
tice to circumvent that as much as possible by imposing 
specifically selected constraints [22, 24 1. Although, even- 
tually, this often leads to rather pleasing results, this 
may not be the case when studying unknown systems 
where good constraints are not known from the outset. 
In addition, since the atomic configuration in RMC is 
not relaxed into a local energy minimum, the resulting 
structure is not necessarily stable. 

The inverse design technique of Franccschctti and 
Zunger allows one to, at least partially, circumvent the 
shortcomings just mentioned by determining the crystal 
structure based on electronic structure properties, which 
are rather sensitive with respect to the atomic positions. 
In their method, an inner local geometry optimization is 
performed in each optimization step to relax the struc- 
ture [25| . However, for the sake of efficiency, the latter 
is conducted using an empirical valence force field only 
[2^. Furthermore, in order to facilitate the calculation, 
they confined themselves to highly symmetric structures 
on a given crystal lattice. 

In this work, we improve upon the existing approaches 
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by proposing a novel and efficient method, which we call 
Inverse Simulated Annealing (ISA). This method com- 
bines the global minimization of a linear combination 
consisting of various geometric and electronic proper- 
ties with structure relaxation to determine an amorphous 
solid in best agreement with available experimental data. 
Specifically, this is achieved by adding the potential en- 
ergy U(R) to the objective function of Eq. [TJ and em- 
ploying a modified hybrid Monte Carlo (HMC)-based SA 
scheme to minimize it. We will demonstrate that the 
present method is efficient enough to be applicable in 
conjunction with accurate electronic structure calcula- 
tions, and in this way allows to routinely determine the 
amorphous structure. 

In the following, we will confine ourselves to effective 
single-particle theories, such as density functional theory 
(DFT) (27| . Hence, the modified objective function to be 
minimized reads as: 

U(R) = U(K) + (X P (R) - X e p xp ) 2 

p 

+ ^ W9 fe[R,{^}]-^) 2 , (2) 



where U (R) is a fictitious and U (R) the potential en- 
ergy, as obtained by DFT, while £ q [R, {ipi}] and ^ xp 
are the computed and experimental values, respectively, 
of an electronic quantity q. 

In minimizing t/(R), we take advantage of the fact that 
by using Eq. ® the accessible phase space is substantially 
reduced and restricted to energetically low-lying configu- 
rations. In other words, even though the dimensionality 
of the phase space is equally vast, the optimization is 
guided in a funnel-like fashion towards the minimum of 
f/(R). Obviously, in spite of that, we still need a global 
optimization method to minimize Eq. [5] that is efficient 
enough to enable the calculation of U(R) at the DFT 
level of theory. The fact that the derivatives of some of 
the properties in Eq. [2] with respect to R are not directly 
available and may not even exist due to possible disconti- 
nuities, immediately suggests a MC-based minimization 
procedure (28|. The development of such a technique is 
therefore an essential part of the present work. 

For the purpose to minimize Eq. ® while at the same 
time using as few as possible electronic structure cal- 
culations, we propose here a novel 'fuzzy' HMC-bascd 
SA scheme within the NVE instead of the more common 
NVT ensemble that consists of only a single modified MD 
step. In comparison to standard MC- or MD-based SA 
techniques, we found that this technique performs partic- 
ularly well as a minimization method, as will be shown. 
The positions and velocities of all atoms in each trial 
move are varied according to a slightly modified velocity- 



Verlet algorithm: 
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where are the ionic velocities, m t the nuclear masses, 
dt a randomly chosen time step from an uniform distribu- 
tion within the interval [0, dt max ], while the prime super- 
scripts are used to indicate quantities of the new (trial) 
configuration. The forces fi are the best possible estimate 
for — dU /dii, i.e. omitting the contributions from those 
terms in the sums of Eq. [2] for which no derivatives are 
directly available. This and the presence of a maximum 
time step dt max , which is in general much larger than 
in standard MD and continuously adjusted to obtain an 
acceptance rate of about 50 %, is why we call our modi- 
fied HMC algorithm 'fuzzy'. In order to ensure that the 
total energy is conserved, in Eq. [3] we have introduced 
an additional prcfactor denoted as C, which chosen in 
such a way that 1/2 £)f mWi? = K ' = E — U' holds, 
where K' is the kinetic energy of the system of the pro- 
posed trial configuration. Within the NVE ensemble, the 
probability of acceptance of a trial move is given by [29[ 



P = min 1, 



U' 



E-U 



37V/2-l^ 



(4) 



where N is the number of atoms. 

As already mentioned, the present approach differs 
from the standard HMC algorithm in the fact that the 
NVE instead of the usual NVT ensemble is employed. 
Furthermore, only a single MD step is taken in each HMC 
step and the velocities are not randomly re-initialized 
thereafter. The necessary random element in our HMC 
method comes from the randomly chosen, variable time 
step dt instead. Whenever a HMC move is accepted, the 
positions and velocities are updated as (r.;, v^) = (r^,v-), 
just as in normal MD. Otherwise, if an HMC move is 
rejected, then one possibility is to maintain (ri,Vj), in 
which case no update is required. We will denote this 
straightforward version of our method as fHMC-NVE. 
However, regarding the efficiency of the minimization 
procedure, it is desirable to design an algorithm that 
combines a large time step with a high acceptance rate. It 
appears that an improvement in this direction is obtained 
by maintaining the velocities of the rejected configura- 
tions, i.e. by updating according to (r^,Vi) = (rj,v^) 
after a rejection. In this modified algorithm, indicated 
hereafter as mfHMC-NVE, the velocities are gradually 
turned in the direction of the forces upon repeated rejec- 
tions. As a consequence, the acceptance probability for 
large displacements (i.e. large dt) increases, since the dis- 
placements become more and more parallel to the forces, 
i.e. the direction of decreasing potential energy. 
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FIG. 1. Comparison of the average over the final potential 
energies at K of amorphous carbon as generated by the 
various minimization method as a function of the quenching 
time At coo i. The averages are based on 40 independent sim- 
ulations, allowing for the calculation of variances and error 
bars as indicated. Note the logarithmic scale for the 'time' 
axis. 



To assess the performance of our HMC-based mini- 
mization technique, we have applied it to carbon using 
the empirical LCBOPII potential [3(|- This bond order 
potential has been shown to accurately describe many 
carbon phases including the disordered, liq uid phase 
within a whole range of different densities 31] . We have 
selected a system consisting of 216 atoms within in a cu- 
bic simulation box with periodic boundary conditions, 
which corresponds to a density of p = 3.1 g/cm 3 , which 
is in close agreement with the experimentally determined 
density of amorphous carbon [33| . For the sake of simplic- 
ity, in these simulations, meant to test and compare the 
performance of different minimization techniques, mo- 
mentarily only the potential energy is minimized. 

The applied total energy schedule as a function of the 
(fictitious) MC 'time', is schematically shown in Fig. QJ,. 
Starting at a high total energy E = —1000 eV, to create 
a well disordered liquid phase, the schedule includes a 
liquid equilibration period at constant E = —1200 eV, 
after which the system is cooled down linearly to E = 
— 1450 eV during a 'time' interval At coo i. After that, the 



system is relaxed in a relatively short quench by further 
decreasing E to a value close to the final potential energy 
Uffi(R). Note that the instantaneous temperature of the 
system can be deduced from K = (3/2)7VfcsT = E — U, 
which implies T = 2(E - U)/(3Nk B ), so that T -> K 
for E -> U ffi . 

The results for the average, final potential energy per 
atom at K, Uffi(R-)/N, as a function of the cooling 
'time' interval At coo i in units of total energy evaluations 
(tee), based on 40 independent simulations, are shown 
in Fig. [TJ; and compared to the results from other, more 
standard minimization techniques. These include the ref- 
erence, a random single atom displacement MC method 
within the NVT ensemble, indicated as RS-MC-NVT, 
and two all atom MC methods within the NVT ensem- 
ble: the RA-MC-NVT method with completely random, 
simultaneous displacements of all atoms and the FA-MC- 
NVT method, where the displacement of each atom is a 
mixture of a random vector and the force on that atom 
with a mixing coefficient chosen such that the efficiency 
is maximized. The applied temperature versus the MC 
'time' for these NVT simulations is schematically given 
in Fig. [TJd. 

As can be seen in Fig. [TJ;, the behavior of the RA- 
MC-NVT technique and the reference is essentially iden- 
tically, which suggests that in the present case it is in- 
significant if either all or a single atom is randomly dis- 
placed. Nevertheless, the straightforward inclusion of nu- 
clear forces in the FA-MC-NVT approach leads to an op- 
timization scheme that can easily get trapped in a local 
minimum and is hence not competitive. On the contrary, 
in the (m)fHMC-NVE method this is circumvented by 
the interplay of dt and C. On average dt max is relatively 
large, i.e. typically about one order of magnitude larger 
than in a conventional MD simulation for carbon and re- 
mains approximately constant during the annealing. In 
this way, the available gradient information is rather well 
exploited. However, upon rejections the decrease of dt is 
counterbalanced by C to conserve the instantaneous total 
energy and therefore prevents the system to be trapped 
in a local minimum. In the end, employing the mfHMC- 
NVE method, the same potential energy than using the 
RS-MC-NVT approach is realized, though with a two 
orders of magnitude shorter cooling time. Comparing 
the final potential energies with the ground state ener- 
gies of diamond (-7.349 eV/atom) and graphite (-7.374 
eV/atom), it is apparent that the eventual structures cor- 
respond to amorphous carbon, whose energies are about 
0.4 eV/atom above the corresponding ground state. 

As already mentioned, the mfHMC-NVE method 
shows the best performance regarding its ability to find 
low energy states. On the other side, it is feasible to 
do much longer simulations (in terms of tee) with the 
RS-MC-NVT method than using the other techniques, 
because the re-evaluation of the total energy after the 
displacement of one single atom is relatively fast for the 
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empirical LCBOPII potential; this is due to the intrin- 
sic local dependencies of the energy contributions in such 
potentials. Since the curve for RA-MC-NVT lies on top 
of that of the RS-MC-NVT technique, the latter is to be 
preferred whenever updating the total energy for single 
atom move is faster than for an all atom move. 
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FIG. 2. Evolution of the (a) potential energy, (b) the sum of 
the squared residuals of the RDF and (c) the tauc-gap as a 
function of time during the optimisation using the mfHMC- 
NVE technique to determine the structure of amorphous car- 
bon. The solid line denotes our novel simulation method, 
while the conventional SA approach is depicted by the dashed 
line. The comparison of the corresponding G(r), as obtained 
using both techniques, with the experimental one [13], given 
by the solid line, is shown in Fig. (2}i. Due to the nearly per- 
fect agreement, the experimental curve is almost completely 
covered by the results from the present method (red line). 

To illustrate our novel method, we apply the mfHMC- 
NVE method to minimize Eq. [2] for amorphous carbon 
at the density functional level of theory (DFT). Therein, 
U(R) is the total energy from DFT supplemented by the 
reduced radial distribution function G(r), derived from 
scattering data, and the optical Tauc gap AE tauc for 
amorphous phases [32|. Hence, in this case, Eq. Stakes 
the form: 

U(R) = U(R) + ™ G £(G„(R) - G^f 



gap (AE Tauc (R) - AE^^f) 2 



(5) 



where G„(R) = G(r n ) denotes a discretized representa- 
tion of G(r), which is defined as G(r) = 47rr(c(r) — c ), 
with c(r) the average (number) density of atoms at a 
distance r and cq the overall density. To obtain a 
smoothened G(r), allowing for the calculation of analyt- 
ical force contributions that were included in the present 
simulations, we have computed it for any r = r„ on a 
grid with a spacing of 0.01 A between the grid points as: 



G(r) 



1 1 

rAr N 



E 



r+ir/2 



-Ar/2 



Pij(r')dr' - Anrco (6) 



where Pij (r) is a Gaussian-shaped polynomial of degree 4 
within the open interval (ry— Ar, r^+Ar) and Pij(r) = 
otherwise, with P^ and dP^jdr being continuous at 
r = rij ± Ar, J Pij(r)dr = 1 and r%j the interatomic 
distance between atom i and j. The values reported for 
the experimental gap of amorphous carbon vary between 
1.0 eV and 2.5 eV, possibly depending on the particular 
sample 18(. Therefore, we have taken an intermediate 
target value equal to AB^f' = 1.7 eV for our simula- 
tion. However, in the present study, we have neglected 
the gradient of the Tauc gap term with respect to R in 
the analytic expression of the forces. Nevertheless, using 
finite differences, it is straightforward to include them, 
although at the price that the computation becomes at 
least a factor of 3N times more expensive. Further de- 
tails on the on-the-fly calculation of the Tauc gap are 
discussed in the Appendix. 

Even though the values of the weight factors wq and 
Wgap have some importance, their impact is relatively 
small. In principle they should be chosen as small as 
possible and just large enough to get a good agreement 
with the experimental data. In the present simulation 
we have used wg = 1 and w gap — 2.5. In general, 
the value Wg should be chosen in such a way that, in 
the beginning of the simulation at high temperature, 
S n (Gn(R) — G>f p ) 2 is on the same order of magnitude 
than the thermal energy ^NksT. In contrast, the pa- 
rameter Wgap can be selected to be considerable smaller 



than 3(AE Tauc (R) - AE^f) 2 /2Nk B T. 

We have linked our code to the CP2K suite of programs 
to compute the necessary total energies and forces [34| . 
The DFT calculations were performed using the Perdew- 
Burke-Ernzcrhof (PBE) exchange correlation functional 



351 ] and norm-conserving Goedecker-type pseudopoten- 
tials [3(|. The total energy schedule applied included an 
equilibration at E = —32950 eV for 1000 tee, followed 
by a cooling from E = -32950 eV to E = -33175 eV 
during 4000 tee and a final run of length 1000 tee during 
which E is further lowered to get as close as possible to 
Uffi- 

The results of such simulations using our novel method, 
with and without the experimental constraints, are pre- 
sented in Fig. [5] and are compared to a conventional 
HMC-bascd SA simulation. The improved overall agree- 
ment with the underlying experimental data is apparent, 
as shown in Figs. and[5Ji. 

We conclude by noting that our novel method in con- 
junction with an appropriate minimization procedure has 
a wide domain of applicability, not limited to amorphous 
phases. We wish to specifically highlight that the present 
scheme can be directly applied to any other disordered 
system, such as liquid water (3?! or, including the NMR 
chemical shift [H, [39| , to determine the structure of pro- 
teins and nucleic acids. Further improvement of the 
method and the minimizer will be presented elsewhere. 
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APPENDIX: TAUC GAP 

The optical Tauc gap AExauc is a convenient definition 
for the gap of amorphous phases, which circumvents the 
difficulty that the band structure of disordered systems 
is not properly defined. It relies on the following relation 
[3^ between the experimental optical gap AE gap and 
the optical absorption coefficient a as a function of the 
photon energy hv: 



a(hv)hv oc (hv — AE gap Y 



(7) 



which is applicable to (amorphous) semi-conductors 
within a certain range of photon energies just beyond 
the gap AEg ap . For the on-the-fly calculation of AExauc 
within each optimization step, we first compute the op- 
tical absorption coefficient a from [4fj| : 



n(E - hv)n(E)dE, 



where K is a constant, while n and n are the densi- 
ties of the occupied and unoccupied states that are com- 
puted from the eigenvalue spectrum of the DFT Hamil- 
tonian after self-consistency has been achieved. Plot- 
ting y Oi(his) hv as a function of hv within a photon en- 
ergy range around the gap obeys a linear regime, from 
which AETauc = AEg a p can be obtained by taking the 
intersection between the linear fit with the horizontal 
y/a(hv) hv = axis. We note that the value of the con- 
stant K is irrelevant for the value of AExauc resulting 
from this approach. Since the linear behavior only ap- 
plies to a finite energy range just beyond the gap, the lin- 
ear fit has to be restricted to this energy interval. For the 
automatic computation of AETauc, we have selected this 
interval to be within the interval ((I — A)fWtot, fWtot), 
where Wtot is the total width of the spectrum. 



* kuehnc@uni-mainz.de 
[1] R. Zallen, The Physics of Amorphous Solids, Wiley, New 
York, 1983. 

[2] S. R. Elliott, Physics of Amorphous Materials, Longman 

Scientific & Technical, Essex, 1990. 
[3] W. H. Zachariasen, J. Am. Chem. Soc. 54, 3841 (1932). 
[4] M. Jansen, J. C. Schon and L. van Wiillen, Angew. 

Chem. Int. Ed. 45, 4244 (2006). 
[5] J. Maddox, Nature 335, 6187 (1988). 
[6] D. M. Deaven and K. M. Ho, Phys. Rev. Lett. 75, 288 

(1995). 

[7] D. J. Wales and H. A. Scheraga, Science 285, 1368 
(1999). 



[8] R. Martonak, A. Laio and M. Parrinello, Phys. Rev. Lett. 

90, 075503 (2003). 
[9] S. Goedecker, J. Chem. Phys. 120, 9911 (2004). 
[10] C. W. Glass, A. G. Oganov and N. Hansen, Comp. Phys. 

Commun. 175, 713 (2006). 
[11] G. Trimarchi and A. Zunger, Phys. Rev. B 75, 104113 

(2007) . 

[12] S. C. Woodley and R. Catlow, Nature Mater. 7, 937 

(2008) . 

[13] C. Wehmeyer, G. F. von Rudorff, S. Wolf, G. Kabbe, D. 
Schaxf, T. D. Kiihne and D. Sebastiani, J. Chem. Phys. 
137, 194110 (2012). 
[14] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, Science 220, 
671 (1983). 

[15] R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 
(1985). 

[16] T. D. Kiihne, M. Krack, F. R. Mohamed and M. Par- 
rinello, Phys. Rev. Lett. 98, 066401 (2007). 
[17] J. Sarnthein, A. Pasquarello and R. Car, Phys. Rev. Lett. 

74, 4682 (2005); Phys. Rev. B 52, 12690 (1995). 
[18] N. A. Marks, D. R. McKenzie, B. A. Pailthorpe, M. 
Bernasconi and M. Parrinello, Phys. Rev. Lett. 76, 768 
(1996); Phys. Rev. B 54, 9703 (1996). 
[19] S. Caravati, M. Bernasconi, T. D. Kiihne, M. Krack and 
M. Parrinello, Appl. Phys. Lett. 91, 171906 (2007). S. 
Caravati, M. Bernasconi, T. D. Kiihne, M. Krack and 
/g\ M. Parrinello, Phys. Rev. Lett. 102, 205502 (2009). S. 

Caravati, M. Bernasconi, T.D. Kiihne, M. Krack and M. 
Parrinello, J. Phys.: Condens. Matter 21, 255501 (2009). 
S. Caravati, D. Colleoni, R. Mazzarello, T. D. Kiihne, 
M. Krack, M. Bernasconi and M. Parrinello, J. Phys.: 
Condens. Matter 23, 265801 (2011). 
[20] M. F. Camellone, T. D. Kiihne and D. Passerone, Phys. 

Rev. B 80, 033203 (2009). 
[21] R. L. McGreevy and L. Pusztai, Mol. Simul. 1, 359 
(1988). 

[22] R. L. McGreevy, J. Phys.: Condens. Matter 13, R877 
(2001). 

[23] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. 

H. Teller and E. Teller, J. Chem. Phys. 21, 1087 (1953). 
[24] S. Kugler, L. Pusztai, L. Rosta, P. Chieux and R. Bellis- 

sent, Phys. Rev B 48, 7685 (1993). 
[25] A. Franceschetti and A. Zunger, Nature 401, 60 (1999). 
[26] P. N. Keating, Phys. Rev. B 149, 674 (1966). 
[27] R. O. Jones and O. Gunnarsson, Rev. Mod. Phys. 61, 
689 (1989). 

[28] S. Duane, A. D. Kennedy, B. J. Pendleton and D. 

Roweth, Phys. Lett. B 2, 216 (1987). 
[29] J. R. Ray, Phys. Rev. A 44, 4061 (1991). 
[30] J. H. Los, L. M. Ghiringhelli, E. J. Meijer and A. Fa- 

solino, Phys. Rev. B 72, 214102 (2005). 
[31] L. M. Ghiringhelli, J. H. Los, A. Fasolino and E. J. Mei- 
jer, Phys. Rev. B 72, 214103 (2005). 
[32] J. Tauc, A. Menth and D. L. Wood, Phys. Rev. Lett. 25, 
749 (1970). 

[33] W. R. Gilkes, P. H. Gaskell, and J. Robertson, Phys. 

Rev. B 51, 12303 (1995). 
[34] J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, 
T. Chassaing and J. Hutter, Comput. Phys. Commun. 
167, 103 (2005). 
[35] J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. 

Lett. 77, 3865 (1996). 
[36] S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 
1703 (1996). 



6 



[37] T. D. Kuhne, M. Krack and M. Parrinello, J. Chem. The- 
ory Comput. 5, 235 (2009). T. D. Kuhne, T. A. Pascal, E. 
Kaxiras and Y. Jung, J. Phys. Chem. Lett. 2, 105 (2011). 
T. A. Pascal, D. Sch&rf, Y. Jung and T. D. Kuhne, J. 
Chem. Phys. 137, 244507 (2012). T. D. Kuhne and R. Z. 
Khaliullin, Nature Commun. 4, 1450 (2013). 

[38] F. Mauri, B. G. Pfrommer and S. G. Louie, Phys. Rev. 
Lett. 77, 5300 (1996). F. Mauri, B. G. Pfrommer and S. 



G. Louie, Phys. Rev. Lett. 79, 2340 (1997). 

[39] D. Sebastiani and M. Parrinello, J. Phys. Chem. A 105, 
1951 (2001). D. Sebastiani and M. Parrinello, Phys. 
Chem. Chem. Phys. 3, 675 (2002). 

[40] R. Bouzerar, C. Amory, A. Zeinert, M. Benlahsen, B. 
Racine, O. Durand-Drouhin and M. Clin, Journal of Non- 
Crystalline Solids 281, 171 (2001). 



